BMC Bioinformatics

External biological knowledge provides valuable information about relationships among metabolites, yet this information is usually not incorporated directly into statistical estimation procedures. Most existing approaches estimate metabolite effects independently, ignoring known biochemical structure such as shared subclasses and pathway membership. We propose a Bayesian hierarchical framework that improves metabolite effect estimates by incorporating external biological information describing relationships among metabolites. The proposed method improves metabolite-specific estimates by allowing related metabolites to borrow information from one another while preserving metabolite-level inference. We evaluate the methodology using simulation studies across a range of sample sizes and heterogeneity regimes together with three metabolomics applications involving distinct biological annotation structures. Across both simulated and real datasets, incorporating external biological information consistently improves metabolite effect estimation. Gains are most pronounced when sample sizes are small and metabolite classes are informative, i.e. more homogenous within classes.

IntroductionCurrent de novo genome assembly tools often demand substantial memory resources, and their execution typically relies on high-performance computing (HPC) clusters. This dependency limits their use in resource-constrained settings. Furthermore, mainstream third-generation sequencing assembly and alignment tools usually require explicit detection of overlap regions between reads, a process that often entails significant computational and storage overhead. ResultsTo address this issue, we developed LightAlign, a lightweight alignment tool for HiFi data that innovatively uses sequence-derived fuzzy features and reduces the peak memory usage during overlap detection. ConclusionsWhen combined with miniasm, LightAlign generated bacterial draft assemblies while maintaining peak memory usage below 1 GB and completed overlap generation for the tested eukaryotic datasets within 1.88 GB RAM.

Motivation: Copy-number variants (CNVs) contribute to human disease and population trait variation. CNV detection from large whole-genome sequencing cohorts remains computationally demanding, as most methods require BAM or CRAM files. Genomic VCF (gVCF) files are smaller, routinely generated by standard variant-calling workflows, and contain the read depth and allelic information needed for CNV detection. However, gVCF files are not directly compatible with established CNV callers that rely on Log R Ratio (LRR) and B Allele Frequency (BAF) signals. Results: We present gVCF2CNV, a Nextflow pipeline that converts gVCF files into Log R Ratio and B Allele Frequency signals compatible with established CNV callers. Applied to 12,509 individuals from the SPARK cohort, gVCF2CNV generated signals at an average of 2.7 million SNV positions per individual and completed signal extraction in 4 hours using 192 CPUs. CNV calling with PennCNV and QuantiSNP identified candidate CNVs across a broad size range, with trio-based Mendelian precision reaching approximately 80% or higher for deletions of at least 30 kb and duplications of at least 5 kb. Application to 414,824 individuals from the All of Us cohort was completed in 96 hours, demonstrating feasibility at biobank scale. These results show that gVCF files can serve as a scalable input for CNV detection in large WGS cohorts. Availability and Implementation: gVCF2CNV is available at https://github.com/JacquemontLab/gVCF2CNV, implemented as a Nextflow pipeline with Perl and Python components, supported on Linux. Contact: mame.seynabou.diop@umontreal.ca

Background: We present Rsearch, an R package that integrates the core functionality of VSEARCH into the R environment and extends it with visualization, parameter optimization, and conversion tools for compatibility with other R packages. By making VSEARCH directly accessible in R, Rsearch lowers the barrier for using VSEARCH and integrating it with downstream statistical and ecological analyses. Results: Comparative analysis with DADA2 using mock community data showed that both pipelines produced relative abundance profiles highly correlated with the expected composition. Compared to DADA2, Rsearch identified fewer OTUs, but these were more consistently prevalent across samples. In contrast, DADA2 appeared to overestimate diversity by splitting sequences into an excessive number of OTUs. In terms of computational performance, vs_cluster_unoise implemented in Rsearch was the fastest of all the clustering and denoising methods, while other Rsearch functions showed runtimes comparable to DADA2. In addition, Rsearch provides functions for systematic optimization of trimming and filtering parameters, an important feature for users who may not otherwise have a clear strategy for parameter selection. The package also includes functions to ensure compatibility with other R packages such as phyloseq. Conclusions: Rsearch offers a practical and accessible framework for analysing metabarcoding data within a single analytical environment and is freely available from The Comprehensive R Archive Network, with the development version hosted on GitHub (https://github.com/CassandraHjo/Rsearch).

Motivation: Real-time qPCR data are commonly analyzed to obtain a Cq values, which are used to compare the relative abundance of a targeted template. However, improper background subtraction distorts Cq and changing the assigned threshold can alter experimental conclusions. In addition, a Cq values provide no information on the quality of the qPCR for a given sample. Alternatively, global fitting of qPCR data provides agnostic quantification without the use of Cq or thresholds, and the fitting terms provide a useful metric of reaction performance. Current implementations of global fitting are limited to custom in-house programs. Results: qPyCR is comprehensive analysis tool written in Python that performs background signal subtractions, computes a threshold, and reports Cq values in a traditional manner. The program also performs global fitting and reports template abundances along with reaction quality values. The quality values can be used to assign performance standards, which allow researchers to distinguish aberrant reactions from true changes in template abundance. Availability and Implementation: Code and example data are freely available on GitHub (https://github.com/sdmoore-labs/qPyCR) and on Zenodo (https://doi.org/10.5281/zenodo.21115342).

Whole genome shotgun sequencing and assembly is routine. However, identifying protein-coding genes in newly assembled genomes remains complex, time-consuming, and labour-intensive. Therefore, most eukaryotic genome assemblies in public databases lack gene annotations reducing their value for evolutionary and functional genomics. Here, we present fDOG-Assembly, a novel tool for targeted, feature architecture-aware ortholog searches directly in unannotated genome assemblies. Benchmarking shows that fDOG-Assembly performs similarly to BUSCO and Compleasm in ortholog identification while offering the advantage of not being restricted to universal single-copy genes. Applied to identify orthologs of 5,000 human genes in rat and Nematostella vectensis, fDOG-Assembly approaches the performance of traditional ortholog search tools that rely on pre-annotated proteomes. Importantly, it can recover orthologs missed by conventional methods because of incomplete gene annotations, helping to fill gaps in phylogenetic profiles. As a case study, we screened 176 soil invertebrate genome assemblies for genes involved in antibacterial compound production. We found that orthologs of {beta}-lactam biosynthesis genes are widespread in springtails, with individual species possessing nearly complete cephamycin biosynthetic gene sets, suggesting they may represent previously unrecognized natural producers of {beta}-lactam antibiotics. Overall, fDOG-Assembly is a powerful resource for orthology-based analyses of the rapidly growing collection of unannotated genome assemblies.

BackgroundSecondary haplotigs retained in a primary assembly inflate genome size and duplicate gene content, and their removal remains a necessary curation step for heterozygous diploid genomes. HapSolo (Solares et al. 2021) treats this step as an optimization problem: it searches three pairwise alignment thresholds (percent identity, query coverage, and the ratio of aligned query to aligned reference length) to minimize a cost defined over single-copy, duplicated, and missing conserved orthologs. The original implementation depended on a separate BUSCO v3 run for every contig, on Python 2.7 and pandas in its inner loop, on a single randomized forward-walk optimizer, and on OrthoDB v9 lineage datasets that are no longer distributed. ResultsHere we report HapSolo2, a reimplementation that replaces the per-contig BUSCO step with miniprot-based ortholog classification against OrthoDB v9, v10, or v12.2 lineage datasets, adds steepest descent and simulated annealing to the original random forward walk, and evaluates all walkers in batched CuPy array operations on a single GPU. On an avocado assembly of 5,122 contigs scored against 1,440 orthologs, 100,000 iterations of random walk and steepest descent ran 5.6- and 7.6-fold faster on one NVIDIA RTX A6000 than on 10 CPU processes. Across 45 runs on four diploid assemblies (Anopheles funestus, Persea americana, Vitis vinifera, and Amblyraja radi-ata; 212 Mb to 3.2 Gb) and three OrthoDB releases, the three optimizers reached best costs within 1.8% of one another. Independent assessment with BUSCO v6.1 showed that HapSolo2 reduced duplicated orthologs from 52.0% to 5.6% in P. americana, from 17.9% to 10.4% in V. vinifera, and from 7.8% to 2.8% in A. funestus, with completeness decreasing by at most 1.3 percentage points. On the avocado assembly, HapSolo2 reduced the cost reached by the original HapSolo from 0.2200 to 0.1036 and raised contig N50 from 3,366 kb to 4,608 kb. ConclusionsHapSolo2 makes ortholog-guided haplotig removal practical in minutes on a single workstation, independent of a BUSCO installation, and reproducible across OrthoDB releases. It is distributed as an installable Python package with a unified command-line pipeline, an assembly statistics and plotting tool, an Apptainer container, and archived copies of the legacy ortholog datasets.

The rapid growth in publicly available genome assemblies has made selecting genomes suitable for downstream analyses increasingly challenging. Differences in assembly and annotation quality can influence gene completeness, duplication rates, contiguity, repeat representation, and other characteristics. Assessing genome quality therefore requires integrating multiple complementary quality metrics that are often generated by independent tools. Here, we present nf-core/genomeqc, a workflow for assessing and comparing genome assemblies. The pipeline accepts RefSeq/GenBank accessions for automatic genome and annotation retrieval, or local genome (FASTA) and annotation (GFF3/GTF) files. It integrates complementary analyses of assembly contiguity, gene completeness, annotation quality, repeat content and other quality metrics using tools such as BUSCO, QUAST, Merqury, and AGAT, before combining the results on a phylogenetic tree for visualisation and comparison across species. GenomeQC is implemented in Nextflow within the nf-core framework, providing an accessible, reproducible, scalable and community-driven workflow for genome quality assessment.

Since copy number variations (CNVs) in pharmacogenes can cause significant alterations in drug metabolism, their reliable detection is of high importance both for large-scale studies and personalized medicine. Whole-genome sequencing, and specifically long-read sequencing, is the gold standard for CNV detection. Despite increasing availability of these technologies, genotyping arrays are still widely used as cost-effective alternatives in biobank and clinical settings, yet calling CNVs based on array intensity signals is challenging due to low base pair resolution. In this work, we developed a neural network model, nnCNV, to predict deletions in the CYP2C19 pharmacogene region from array intensity signals. We compared our method to the most widely used algorithm, PennCNV, and demonstrated better performance reaching 100% accuracy in the test dataset. Furthermore, we predicted probe-by-probe CYP2C19 deletion coordinates for all Estonian Biobank samples using nnCNV and PennCNV, and validated these predictions using an identity-by-descent (IBD) sharing method, which also demonstrated superior nnCNV performance. For the deletion samples with conflicting PennCNV and nnCNV predictions, we performed PCR analysis for validation, which showed 97% precision for nnCNV compared to 23% for PennCNV. Finally, we assessed the gradient-based feature importance maps and showed that nnCNV utilizes signal intensity information not only from deletion probes, but also from probes in flanking regions. Our results demonstrate that long-range information, which cannot be utilized by hidden Markov models, can improve CNV calling.

MotivationHigh-dimensional microarray datasets remain valuable for cancer biomarker discovery, but their small sample sizes make robust and interpretable feature selection challenging. Efficient workflows are needed to derive compact gene signatures while preserving biological interpretability. ResultsWe developed a two-stage biomarker-discovery workflow that combines cross-validated XGBoost rank aggregation with support vector machine recursive feature elimination and cross-validation (SVM-RFECV) to identify compact candidate biomarker panels. The workflow was evaluated on 21 public binary and multiclass microarray datasets using repeated stratified cross-validation for internal validation. Across the dataset collection, the selected panels demonstrated strong internal discriminative performance while remaining sufficiently compact for downstream biological interpretation. SHAP analysis identified dataset- and class-specific discriminative genes, and functional enrichment analysis supported the biological coherence of representative consensus signatures. The proposed workflow provides an interpretable and reproducible framework for candidate biomarker discovery from small-sample microarray datasets. AvailabilitySource code and processed outputs are freely available at https://github.com/mashiyat-mahjabin-prapty/microarray-feature-selection.

Background. Cleavage under targets and release using nuclease (CUT&RUN) maps the genome-wide locations of chromatin-associated proteins and provides an improved alternative to chromatin immunoprecipitation sequencing (ChIP-seq) for profiling sequence-specific transcription factor binding sites. Identifying these binding sites plays a critical role in understanding gene regulation, and transcription factors provide a useful setting for benchmarking because their well-defined sequence motifs serve as built-in controls for evaluating performance. Compared with ChIP-seq, CUT&RUN achieves higher resolution and lower background by avoiding cross-linking and bulk precipitation. Its distinct fragment length and cleavage characteristics, however, limit the direct transfer of existing computational tools, which primarily target ChIP-seq data. The performance of these tools on CUT&RUN can depend strongly on preprocessing choices. In this work, we investigate preprocessing strategies for transcription factor CUT&RUN, focusing on fragment length filtering and spike-in calibration. We aim to improve peak detection and provide practical guidance for analysis. Results. We designed a benchmarking method to evaluate peak-calling procedures for CUT&RUN data and the effects of preprocessing approaches, including fragment length filtering and spike-in calibration. We benchmarked the two most widely used peak callers, MACS2 and SEACR, by assessing motif enrichment---the degree to which identified peaks contain the expected transcription factor binding motifs. Filtering for fragments with a length [≤]120 bp generally improved target motif enrichment. Spike-in calibration using heterologous Saccharomyces cerevisiae DNA improved motif elucidation substantially for MACS2, with little benefit for SEACR. By contrast, using Escherichia coli DNA as a spike-in control often failed to produce valid results unless we could meticulously control E. coli contamination. MACS2 performed robustly across samples. SEACR performed especially well on clean, sparse-background datasets, but performed poorly on some datasets with denser background signal and often produced numerous apparent false positives. While MACS2 provided robust results under minor perturbations in fragment length filtering, SEACR exhibited greater sensitivity to such changes. Discussion. Our benchmarking highlights how both peak caller choice and preprocessing strategy shape the analysis of transcription factor CUT&RUN data. By comparing the robustness and limitations of two widely used peak callers, we provide practical guidance on fragment length filtering, spike-in calibration, and tool selection. These findings help improve the processing and interpretation of CUT&RUN data, allowing researchers to more rapidly and reliably utilize this new technology. We expect that our work will guide more informed choices in CUT&RUN analysis and support the development of improved computational methodologies.

Gene Ontology (GO) Biological Process overrepresentation analysis is widely used to interpret gene lists from genetic studies, yet results depend critically on the background (universe/reference list) against which enrichment is tested. This paper examines how genome-exome background mismatch alters GO Biological Process significance and induces annotation-driven bias. First, Monte Carlo simulations across multiple input gene list sizes show that enrichment p-values shift systematically when lists sampled from an exome-like universe are tested against a genome background (and vice versa), producing both inflation and deflation of significance depending on GO term composition; these shifts increase with gene list size. Second, applied analyses of gene lists derived from Genome-Wide Association Studies (GWAS) and Whole Exome Studies (WES) across brain, immune, and metabolic domains demonstrate that background choice changes the set of significant GO IDs, yielding reference-specific terms consistent with both Type I errors (false positives) and Type II errors (false negatives). Because genome backgrounds are commonly used by default, the practical risk is greatest when WES-derived lists are analyzed with genome reference lists. To support reproducible best practice, we provide a simple command set for selecting and documenting study-appropriate backgrounds and for assessing sensitivity of GO Biological Process results to the chosen universe.

We have developed a method to integrate transcriptomic coexpression networks across diverse experimental conditions within a single species. Our framework expands previous pannetwork approaches into a hypergraph based pannetwork. It first identifies coexpressed gene clusters within each individual dataset and extracts frequently co-expressed gene sets across multiple datasets using a frequent itemset mining algorithm. This process yields a hypergraph where each hyperedge is assigned a frequency (universality, U). We then extracted a subnetwork comprising high U hyperedges as the core network and formed modules within it. We applied our method to 106 Escherichia coli transcriptome datasets from the GEO database. The modularity of the core network peaked at a universality cutoff of 15, which was subsequently used to define it. Approximately 70% of the resulting core modules correspond to operons, and conversely, approximately 70% of all operons are covered by these core modules. We visualized the core network via an inter-modular network and analyzed core module dataset relationships using a modularity profile matrix. Based on these analyses, we successfully visualized the dynamic reorganization of the coexpression network in response to environmental changes in bacteria.

In recent years, the issue of reproducibility in scientific studies has been at the forefront of scientific discussion. Many of these issues arise due to poor employment of statistical methods, making it crucial to address these knowledge gaps for frequently used methods. Elastic net regression (ENET) is one such algorithm used frequently in biology and public health due to its ability to select important predictors (also known as "features") from large data sets. However, ENET struggles with consistent feature selection and hyperparameter choice in these applications due to low observation counts and noisy data. This limitation hinders ENET's effectiveness in biological research due to low model reproducibility and predictor interpretability. We propose the Instability Quotient Criterion (IQC) as a metric to quantify the consistency of selected features for feature selection models such as ENET, with a specific emphasis on biological applications. IQC leverages ensemble modeling and Shannon's entropy in order to determine the stability of selected features derived from ENET models. We demonstrate IQC's robustness via simulation studies, showing that IQC is capable of detecting significant decreases in model stability with lower observation counts, higher dimensionality, and higher proportions of noisy signal. We additionally utilize simulation to provide a naive threshold for low versus high stability. Finally, we applied IQC to a case study to illustrate how our method can be integrated into a larger analysis pipeline when looking at microarray gene expression data. IQC indicated stability issues for differentially-expressed genes when classifying acute leukemia subtypes in spite of high classification success. Lastly, we discuss how model instability impedes interpretation of biologically significant genes for discerning cancer subtypes. Overall, we present IQC as a robust evaluation metric that can easily be incorporated into high-dimensional biological analyses.

Genetic association analyses of exome sequencing data pose distinct methodological challenges in small cohorts, which are common in rare disease research. Data filtering prior to association testing requires particular care, especially when analyzing datasets generated under different study designs. Most of current genetic association tools need to be integrated into complex data processing pipelines, whose implementation requires technical expertise and is prone to error. To alleviate mentioned issues we built a comprehensive, high performance Nextflow workflow for genetic association testing, which includes all commonly performed data filtering and preparation steps, controlled by a rich set of modifiable parameters. We compared our software against existing tools incorporating quality control and aggregated association testing steps and found our solution to offer the most complete suite of data preparation, association testing and reporting modules, making small-cohort association analyses accessible to teams with limited technical expertise.

Rapid and accurate microbial identification is critical for interpreting biological data in basic research and when making applied decisions on how to effectively treat patients and control human, animal, and plant diseases. Advancements in high-throughput sequencing have the potential to expedite fungal species identification and thus fungal biological research; however, analyses of ever-increasing numbers of genomes also present computational challenges. In this work, we evaluated how whole-genome similarity can serve as the basis for accurate classification and identification across the Kingdom Fungi and the Phylum Oomycota. Results from the computationally efficient k-mer-based tool sourmash are compared with those from more computationally demanding BLAST-based similarity method ANIb, as well as with conventional phylogenomic approaches, including maximum-likelihood concatenated ortholog trees and SNP-based methods. We observed that sourmash delivers orders-of-magnitude gains in speed and memory efficiency while maintaining strong concordance with phylogenomic methods. We found that a k-mer size of 21 is robust for genus and species determination, and that larger k-mers are well-suited for identification below the species rank. These results demonstrate that k-mer-based whole-genome similarity provides a scalable framework for fungal classification, enabling rapid analysis, lowering computational and bioinformatic barriers, and supporting the development of efficient identification pipelines.

Comparisons of genome structure between related organisms are central to understanding genome evolution, gene family dynamics, and the genomic basis of phenotypic variation. Synteny, the conserved co-localization of genes along chromosomes, is most readily interpreted visually, yet many existing synteny visualization tools require local software installation, command-line proficiency, and/or dedicated server infrastructure, and produce static images that cannot be explored interactively. Here, I present SyntenyPair Explorer, a lightweight, installation-free tool for interactive visualization of synteny between two genomes. The application runs entirely within a standard web browser as a single, self-contained HTML file with no external dependencies and no server-side component. SyntenyPair Explorer accepts standard file formats already produced by common comparative genomics workflows, including FASTA genome assemblies (used to compute optional assembly summary statistics), GFF3/GTF gene annotations, and either BLAST tabular output (outfmt 6) or MCScanX collinearity files. Syntenic relationships are resolved by gene-identifier matching between the relationship file and the gene annotations, so that each relationship corresponds to a discrete gene-to-gene link. Interactive features include continuous zoom and pan, gene search with automatic centering of the partner genome on the syntenic counterpart, synteny block coloring, extensively customizable gene highlights and annotation callouts, session saving and restoration, and publication-quality image export. I demonstrate the tool by visualizing structural differences at the alpha-amylase loci between two strains of the industrially important fungus Aspergillus oryzae. SyntenyPair Explorer lowers the technical barrier to interactive synteny visualization and is freely available under the MIT license at https://github.com/GibbonsLabGenomics/SyntenyPair-Explorer, with a live browser-based version at https://gibbonslabgenomics.github.io/SyntenyPair-Explorer/.

1 Abstract 1.1 Background Many longitudinal omics studies contain only a small number of repeated measurements collected before, during, or after an intervention. Existing approaches, including mixed-effects models and generalized additive models, estimate temporal effects but do not generally provide a discrete representation of trajectory topology that can be queried directly across experimental groups. 1.2 Methods We developed LongOmicsTraj, an open-source R package for topology-based representation and querying of short longitudinal omics trajectories. The framework encodes the direction of change between adjacent visits as up, down, or flat, with the ordered sequence defining an Ordinal Trajectory State (OTS). LongOmicsTraj operates downstream of trajectory estimation and can therefore be applied to empirical summaries or model-derived visit-level estimates, including those from linear mixed-effects models, generalized additive models, and polynomial regression, following a maSigPro-style time-course formulation [1]. OTS labels provide a common representation for topology-based querying, cross-group comparison, and evaluation of higherlevel representations such as trajectory clusters. We evaluated the framework using controlled simulations and bronchial biopsy transcriptomic data from the GLUCOLD corticosteroid intervention study (GEO accession GSE36221), measured at baseline, 6 months, and 30 months. The biological analysis compared continued inhaled corticosteroid (ICS) treatment, ICS withdrawal after 6 months, and placebo. 1.3 Results In simulations, LongOmicsTraj recovered predefined stable, monotonic, transient, rebound, and oscillatory trajectories with high accuracy when longitudinal signal was sufficiently clear, with performance declining under high-noise conditions and depending partly on the upstream estimator. In GLUCOLD, comparator-aware topology queries reduced 20,358 measured transcripts to 168 genes showing a corticosteroid response that was maintained during continued treatment, reversed following withdrawal, and was not reproduced under placebo. The selected genes included established corticosteroid-response genes and were enriched for immune-cell migration, chemotaxis, cell adhesion, and extracellular-matrix organisation. Topology-aware evaluation of FlexMix trajectory clusters additionally revealed substantial within-cluster temporal heterogeneity, with topology purities of approximately 46% to 60%. 1.4 Conclusions LongOmicsTraj provides a compact, directly queryable representation of temporal direction and order in short longitudinal omics studies. It complements existing longitudinal estimation and clustering methods by making trajectory structure explicit, enabling structured cross-group queries and quantification of temporal heterogeneity within trajectory clusters.

High-throughput single-cell sequencing provides a scalable solution for characterizing cells and profiling gene expression for hundreds to millions of cells. However, this process gives rise to doublets, which can lead to inaccurate conclusions drawn from the data. A number of packages have therefore been developed to help accurately detect them, and in particular scDblFinder has been shown to outperform alternatives in the detection of doublets in single-cell (RNA) sequencing data. Being implemented in R, however, its adoption has been more limited in the Python community. Here, we present scDblFinderPy, a Python-based implementation of the scDblFinder R method, and show that it obtains similar performances. Furthermore, we include in it optional GPU support, thus further speeding up the process.

Whole plate detection of calcium flux signals using fluorescence is a powerful technology for monitoring oscillatory activity of excitable cells in a high throughput fashion that sees broad applications in neurotoxicity and cardiotoxicity research. However, analysis of calcium oscillation profiles can output more than 60 peak kinetic parameters, which makes manual data analysis challenging. We have created the CalFluxTools R package, an automated pipeline for parsing plate data for 384 well plates, generating quality control metrics and plots, and performing data analysis, including T-tests, Z-score analysis, and machine learning prediction of compound toxicity values. CalFluxTools requires a plate map and experimental manifest that includes parameter filtration specifications, as well as data in the ".statall" format. CalFluxTools runs in less than a minute on a Macbook Pro with a 2.6 GHz 6-Core Intel Core i7 processor.